function res = residual_for_ztilde(Par,phi,vz,agg,vlabor,logztildeH)
% Description: 
%   Returns the residuals of the first order condition of the union's 
%   problem for given levels of the cutoff productivity, taking labor 
%   allocation and aggregate variables as given
% Output:
%   res - scaler for residual
% Inputs: 
%   Par - parameter structure
%   phi - scaler, parameter phi for the Home firm
%   vz - 2 by 1 vector of ex ante productivity ([zHome;zForeign])
%   agg - 1 by 3 vector of equilibrium values of aggregate variables
%   vlabor - 4 by 1 vector of labor allocation
%   logztildeH - scaler, guess of the cutoff productivity
% Note: Cutoff productivity in log form (logztildeH) to accommodate 
%   searches of negative numbers by fsolve
% Note: Productivity and labor input CANNOT be matrices; Only vectors (for
%   individual sectors) are allowed for this function

    % Load and rename variables
    sigma = Par.sigma;
    rho = Par.rho;
    tau = Par.tau;
    kappa = Par.kappa;
    sigmazH = Par.sigmazH;
    zbarH = log(vz(1,1))-sigmazH^2/2;
    zF = vz(2,1);
    ztildeH = exp(logztildeH);
    % Compute origin-destination level output, prices, and profits by strike outcomes
    y_nostrike = vlabor.*kron([ztildeH;zF],[1;1])./[1;tau;tau;1];
    y_strike = y_nostrike.*kron([1-kappa;1],[1;1]);
    ysector_nostrike = [(y_nostrike(1,:).^((rho-1)/rho)+y_nostrike(3,:).^((rho-1)/rho)).^(rho/(rho-1));(y_nostrike(2,:).^((rho-1)/rho)+y_nostrike(4,:).^((rho-1)/rho)).^(rho/(rho-1))];
    ysector_strike = [(y_strike(1,:).^((rho-1)/rho)+y_strike(3,:).^((rho-1)/rho)).^(rho/(rho-1));(y_strike(2,:).^((rho-1)/rho)+y_strike(4,:).^((rho-1)/rho)).^(rho/(rho-1))];
    psector_nostrike = ([agg(1);agg(2)]./ysector_nostrike).^(1/sigma);
    psector_strike = ([agg(1);agg(2)]./ysector_strike).^(1/sigma);
    p_nostrike = repmat(psector_nostrike.*ysector_nostrike.^(1/rho),2,1)./y_nostrike.^(1/rho)./[1;tau;tau;1];
    p_strike = repmat(psector_strike.*ysector_strike.^(1/rho),2,1)./y_strike.^(1/rho)./[1;tau;tau;1];
    profits_nostrike = sum([1;tau].*p_nostrike(1:2,1).*y_nostrike(1:2,1) - vlabor(1:2,1),1);
    profits_strike = sum([1;tau].*p_strike(1:2,1).*y_strike(1:2,1) - (1-kappa).*vlabor(1:2,1),1);
    % Compute the marginal revenue product of productivity by whether strikes take place
    mrpz_nostrike = fun_mrpz(Par,agg,[ztildeH;zF],vlabor);
    mrpz_strike = fun_mrpz(Par,agg,[(1-kappa)*ztildeH;zF],vlabor);
    % Calculate the D function as defined by (47) in the notes
    Dz = 1/(sum(mrpz_nostrike(1:2,1),1) - (1-phi)*(1-kappa)*sum(mrpz_strike(1:2,1),1));
    % Compute the right-hand side of (49) in the notes
    RHS = Dz*(kappa*(vlabor(1)+vlabor(2))+profits_nostrike-profits_strike);
    % Compute the left-hand side of (49) in the notes
    LHS = (1-logncdf(ztildeH,zbarH,sigmazH))/lognpdf(ztildeH,zbarH,sigmazH);
    % Calculate the percentage deviation of the RHS from the LHS of (49)
    res = RHS/LHS - 1;
end